Лабораторная работа 2¶
Сыромятников Дмитрий (КН-301)
Решить задачу Коши: $$ \begin{cases} \frac{dy}{dx} = 30y(x - 0.2)(x - 0.7) \\ y(0) = 0.1 \\ x \in [0, 1] \end{cases} $$
используя различные методы решения диффференциальных уравнений:
- Явный метод Эйлера
- Метод Эйлера с пересчетом
- Метод Коши
- Метод Рунге-Кутта 4-го порядка
- Неявный метод Эйлера
- Метод Тейлора 2 порядка точности на одном шаге
- Метод Тейлора 3 порядка точности на одном шаге
- 2-хшаговый метод Адамса
- 3-хшаговый метод Адамса
Решения для каждого метода построить по 5, 10, 50 и 100 узлам.
Для применения методов Адамса осуществлять разгон методом Рунге-Кутта 4-го порядка.
Вычислить точное решение задачи. На графике отобразить ломанные, построенные перечисленными методами, точное решение. На основе графиков сделать выводы о погрешностях вычислений у перечисленных методов.
Решение¶
Вычисление точного решения¶
Общее решение
$\frac{dy}{dx} = 30 \cdot (x - 0.2)(x - 0.7)$ $| \cdot \frac{dx}{dy}$
$\frac{dy}{y} = 30(x - 0.2)(x - 0.7)dx$ $| \int$
$\int \frac{dy}{y} = \int 30(x - 0.2)(x - 0.7)dx$
$\ln |y| = 30(\frac{x^3}{3} - 0.9 \cdot \frac{x^2}{2} + 0.14 \cdot {x}) + C$
$\ln |y| = 10x^3 - 13.5x^2 + 4.2x + C$
$\Rightarrow y = C \cdot \exp (10x^3 - 13.5x^2 + 4.2x)$
Частное (искомое) решение
$y(0) = C \cdot 1 = 0.1 \rightarrow C = 0.1$
$\Rightarrow y(x) = 0.1 \cdot \exp (10x^3 - 13.5x^2 + 4.2x)$
from math import exp
def y_exact(x):
"""
Точное решение дифференциального уравнения
:param x: точка
:return: значение решения в точке
"""
return 0.1 * exp(10 * x ** 3 - 13.5 * x ** 2 + 4.2 * x)
Для вычисления решения методами Тейлора требуются частные производные правой части данного дифференциального уравнения ($f^{'}_x, f^{'}_y$):
$f(x, y) = 3y(x - 0.2)(x - 0.7)$
$f^{'}_x(x, y) = 6xy - 2.7y$
$f^{'}_y(x, y) = 3(x - 0.2)(x - 0.7)$
def f(x, y):
"""
Функция f данного дифференциального уравнения
"""
return 30 * y * (x - 0.2) * (x - 0.7)
def df_dx(x, y):
"""
Частная производная функции f данного дифференциального уравнения по x
"""
return 6 * x * y - 2.7 * y
def df_dy(x, y):
"""
Частная производная функции f данного дифференциального уравнения по y
"""
return 3 * (x - 0.2) * (x - 0.7)
Методы численного решения ОДУ¶
def _get_h(a, b, m):
"""
:param a: начало отрезка
:param b: конец отрезка
:param m: число узлов по которым следует разбить отрезок
:return: расстояние между соседними узлами
"""
return (b - a) / (m - 1)
def euler_explicit(f, a, b, y0, n):
"""
Решает задачу Коши для обыкновенного дифференциального уравнения (ОДУ)
первого порядка явным методом Эйлера
:param f: Функция правой части ОДУ вида dy/dx = f(x, y)
:param a: Начальная точка интервала интегрирования
:param b: Конечная точка интервала интегрирования (должна быть больше a)
:param y0: Начальное условие (значение y в точке a)
:param n: Количество узлов сетки на отрезке [a, b]
:return: Кортеж списков X, Y. X - узлы, Y - вычисленные значения функции в узлах
"""
h = _get_h(a, b, n)
X = [a + i * h for i in range(n)]
Y = [y0]
for i in range(n - 1):
Y.append(Y[i] + h * f(X[i], Y[i]))
return X, Y
def euler_recalc(f, a, b, y0, n):
"""
Решает задачу Коши для обыкновенного дифференциального уравнения (ОДУ)
первого порядка методом Эйлера с пересчетом
:param f: Функция правой части ОДУ вида dy/dx = f(x, y)
:param a: Начальная точка интервала интегрирования
:param b: Конечная точка интервала интегрирования (должна быть больше a)
:param y0: Начальное условие (значение y в точке a)
:param n: Количество узлов сетки на отрезке [a, b]
:return: Кортеж списков X, Y. X - узлы, Y - вычисленные значения функции в узлах
"""
X, Yt = euler_explicit(f, a, b, y0, n)
h2 = _get_h(a, b, n) / 2
Y = [y0]
for i in range(1, n):
Y.append(Y[i - 1] + h2 * (f(X[i - 1], Y[i - 1]) + f(X[i], Yt[i])))
return X, Y
def cauchy(f, a, b, y0, n):
"""
Решает задачу Коши для обыкновенного дифференциального уравнения (ОДУ)
первого порядка явным методом Коши
:param f: Функция правой части ОДУ вида dy/dx = f(x, y)
:param a: Начальная точка интервала интегрирования
:param b: Конечная точка интервала интегрирования (должна быть больше a)
:param y0: Начальное условие (значение y в точке a)
:param n: Количество узлов сетки на отрезке [a, b]
:return: Кортеж списков X, Y. X - узлы, Y - вычисленные значения функции в узлах
"""
X2, Y2 = euler_explicit(f, a, b, y0, 2 * n - 1)
h = _get_h(a, b, n)
Y = [y0]
for i in range(1, n):
k2 = 2 * i - 1
Y.append(Y[i - 1] + h * f(X2[k2], Y2[k2]))
return X2[::2], Y
def runge_kutta(f, a, b, y0, n):
"""
Решает задачу Коши для обыкновенного дифференциального уравнения (ОДУ)
первого порядка методом Ренге-Кутта 4-го порядка
:param f: Функция правой части ОДУ вида dy/dx = f(x, y)
:param a: Начальная точка интервала интегрирования
:param b: Конечная точка интервала интегрирования (должна быть больше a)
:param y0: Начальное условие (значение y в точке a)
:param n: Количество узлов сетки на отрезке [a, b]
:return: Кортеж списков X, Y. X - узлы, Y - вычисленные значения функции в узлах
"""
h = _get_h(a, b, n)
X = [a + h * i for i in range(n)]
Y = [y0]
for i in range(1, n):
Xi = X[i - 1]
Yi = Y[i - 1]
K1 = f(Xi, Yi)
K2 = f(Xi + h / 2, Yi + h * K1 / 2)
K3 = f(Xi + h / 2, Yi + h * K2 / 2)
K4 = f(Xi + h, Yi + h * K3)
Y.append(Yi + (h / 6) * (K1 + 2 * K2 + 2 * K3 + K4))
return X, Y
def euler_implicit(f, a, b, y0, n):
"""
Решает задачу Коши для обыкновенного дифференциального уравнения (ОДУ)
первого порядка неявным методом Эйлера
:param f: Функция правой части ОДУ вида dy/dx = f(x, y)
:param a: Начальная точка интервала интегрирования
:param b: Конечная точка интервала интегрирования (должна быть больше a)
:param y0: Начальное условие (значение y в точке a)
:param n: Количество узлов сетки на отрезке [a, b]
:return: Кортеж списков X, Y. X - узлы, Y - вычисленные значения функции в узлах
"""
X, Ye = euler_explicit(f, a, b, y0, n)
Y = [y0]
h = _get_h(a, b, n)
for i in range(1, n):
Y.append(Y[i - 1] + h * f(X[i], Ye[i]))
return X, Y
def taylor2(f, a, b, y0, n):
"""
Решает задачу Коши для обыкновенного дифференциального уравнения (ОДУ)
первого порядка методом Тейлора 2-го порядка
:param f: Функция правой части ОДУ вида dy/dx = f(x, y)
:param dfdx: Частная производная f по x
:param dfdy: Частная производная f по y
:param d2fdxdx: Вторая частная производная f по x
:param d2fdydy: Вторая частная производная f по y
:param d2fdxdy: Вторая смешанная производная f
:param a: Начальная точка интервала интегрирования
:param b: Конечная точка интервала интегрирования (должна быть больше a)
:param y0: Начальное условие (значение y в точке a)
:param n: Количество узлов сетки на отрезке [a, b]
:return: Кортеж списков X, Y. X - узлы, Y - вычисленные значения функции в узлах
"""
return euler_explicit(f, a, b, y0, n)
def taylor3(f, dfdx, dfdy, a, b, y0, n):
"""
Решает задачу Коши для обыкновенного дифференциального уравнения (ОДУ)
первого порядка методом Тейлора 3-го порядка
:param f: Функция правой части ОДУ вида dy/dx = f(x, y)
:param dfdx: частная производная f по x
:param dfdy: частная производная f по y
:param a: Начальная точка интервала интегрирования
:param b: Конечная точка интервала интегрирования (должна быть больше a)
:param y0: Начальное условие (значение y в точке a)
:param n: Количество узлов сетки на отрезке [a, b]
:return: Кортеж списков X, Y. X - узлы, Y - вычисленные значения функции в узлах
"""
h = _get_h(a, b, n)
X = [a + i * h for i in range(n)]
Y = [y0]
for i in range(0, n - 1):
Xi = X[i]
Yi = Y[i]
y = Yi + h * f(Xi, Yi) + (h * h / 2) * (dfdx(Xi, Yi) + dfdy(Xi, Yi) * f(Xi, Yi))
Y.append(y)
return X, Y
def adams2(f, a, b, y0, n):
"""
Решает задачу Коши для обыкновенного дифференциального уравнения (ОДУ)
первого порядка методом 2-шаговым методом Адамса. Разгон производит, используя метод
Рунге-Кутта 4-го порядка
:param f: Функция правой части ОДУ вида dy/dx = f(x, y)
:param a: Начальная точка интервала интегрирования
:param b: Конечная точка интервала интегрирования (должна быть больше a)
:param y0: Начальное условие (значение y в точке a)
:param n: Количество узлов сетки на отрезке [a, b]
:return: Кортеж списков X, Y. X - узлы, Y - вычисленные значения функции в узлах
"""
h = _get_h(a, b, n)
X = [a + i * h for i in range(n)]
Y = [y0]
X0 = X[0]
Y0 = Y[0]
K1 = f(X0, Y0)
K2 = f(X0 + h / 2, Y0 + h * K1 / 2)
K3 = f(X0 + h / 2, Y0 + h * K2 / 2)
K4 = f(X0 + h, Y0 + h * K3)
Y.append(Y0 + (h / 6) * (K1 + 2 * K2 + 2 * K3 + K4))
for i in range(2, n):
yi = Y[i - 1] + (h / 2) * (3 * f(X[i - 1], Y[i - 1]) - f(X[i - 2], Y[i - 2]))
Y.append(yi)
return X, Y
def adams3(f, a, b, y0, n):
"""
Решает задачу Коши для обыкновенного дифференциального уравнения (ОДУ)
первого порядка методом 3-шаговым методом Адамса. Разгон производит, используя метод
Рунге-Кутта 4-го порядка
:param f: Функция правой части ОДУ вида dy/dx = f(x, y)
:param a: Начальная точка интервала интегрирования
:param b: Конечная точка интервала интегрирования (должна быть больше a)
:param y0: Начальное условие (значение y в точке a)
:param n: Количество узлов сетки на отрезке [a, b]
:return: Кортеж списков X, Y. X - узлы, Y - вычисленные значения функции в узлах
"""
h = _get_h(a, b, n)
X = [a + i * h for i in range(n)]
Y = [y0]
for i in range(2):
Xi = X[i]
Yi = Y[i]
k1 = f(Xi, Yi)
k2 = f(Xi + h/2, Yi + h * k1 / 2)
k3 = f(Xi + h/2, Yi + h * k2 / 2)
k4 = f(Xi + h, Yi + h * k3)
y_next = Yi + (h / 6) * (k1 + 2*k2 + 2*k3 + k4)
Y.append(y_next)
for i in range(3, n):
yi = Y[i - 1] + (h / 12) * (23 * f(X[i - 1], Y[i - 1]) - 16 * f(X[i - 2], Y[i - 2]) + 5 * f(X[i - 3], Y[i - 3]))
Y.append(yi)
return X, Y
Параметры¶
# Параметры задачи
a, b = 0, 1
y01 = 0.1
ns = [5, 10, 50, 100]
# ns = range(5, 100, 2) # важно начинать с 5 узлов
hs = [_get_h(a, b, n) for n in ns]
Данные¶
# Точные значения функции
X_exact = [0 + 0.001*i for i in range(1001)]
Y_exact = [y_exact(x) for x in X_exact]
# Расчет всех методов для всех n
data = {}
for n in ns:
Xe, Ye = euler_explicit(f, a, b, y01, n)
Xer, Yer = euler_recalc(f, a, b, y01, n)
Xc, Yc = cauchy(f, a, b, y01, n)
Xrc, Yrc = runge_kutta(f, a, b, y01, n)
Xei, Yei = euler_implicit(f, a, b, y01, n)
Xt2, Yt2 = taylor2(f, a, b, y01, n)
Xt3, Yt3 = taylor3(f, df_dx, df_dy, a, b, y01, n)
Xa2, Ya2 = adams2(f, a, b, y01, n)
Xa3, Ya3 = adams3(f, a, b, y01, n)
data[n] = {
'euler_explicit': (Xe, Ye),
'euler_recalc': (Xer, Yer),
'cauchy': (Xc, Yc),
'runge_kutta': (Xrc, Yrc),
'euler_implicit': (Xei, Yei),
'taylor2': (Xt2, Yt2),
'taylor3': (Xt3, Yt3),
'adams2': (Xa2, Ya2),
'adams3': (Xa3, Ya3)
}
import plotly.graph_objects as go
from plotly.subplots import make_subplots
fig = make_subplots()
fig.add_trace(go.Scatter(
x=X_exact, y=Y_exact,
name='Точное решение',
line=dict(width=4, color='#000000'),
hovertemplate='%{y:.4f}<extra></extra>'
))
methods = [
('euler_explicit', 'Явный метод Эйлера', '#E6194B'), # Ярко-красный
('euler_recalc', 'Метод Эйлера с пересчетом', '#3CB44B'), # Зеленый "кричащий"
('cauchy', 'Метод Коши', '#FFE119'), # Ярко-желтый
('runge_kutta', 'Метод Рунге-Кутта', '#4363D8'), # Интенсивный синий
('euler_implicit', 'Неявный метод Эйлер', '#F58231'), # Оранжевый
('taylor2', 'Метод Тейлора 2 порядка', '#00FF7F'),
('taylor3', 'Метод Тейлора 3 порядка', '#FFC0CB'),
('adams2', '2-шаговый метод Адамса', '#911EB4'), # Фиолетовый
('adams3', '3-шаговый метод Адамса', '#42D4F4') # Бирюзовый (циан)
]
for method_key, name, color in methods:
X, Y = data[5][method_key]
fig.add_trace(go.Scatter(
x=X, y=Y,
name=name,
line=dict(width=2, color=color),
hovertemplate='%{y:.4f}<extra></extra>'
))
frames = []
for n in ns:
frame_data = []
for method_key, name, color in methods:
X, Y = data[n][method_key]
frame_data.append(
go.Scatter(x=X, y=Y, line=dict(color=color))
)
frames.append(go.Frame(
data=frame_data,
name=str(n),
traces=[i for i in range(1, len(data[5]) + 1)]
))
fig.frames = frames
sliders = [{
'active': 0,
'currentvalue': {'prefix': 'Число узлов (n): '},
'steps': [
{
'args': [
[str(n)],
{
'frame': {'duration': 0, 'redraw': True},
'mode': 'immediate',
'title.text': f'Сравнение методов (n={n}, h={h_val:.4f})'
}
],
'label': f'n={n} (h={h_val:.4f})',
'method': 'animate'
} for n, h_val in zip(ns, hs)
]
}]
# Обновление макета с динамическим заголовком
fig.update_layout(
title=f'Сравнение методов решения ОДУ (n={ns[0]}, h={hs[0]:.4f})',
xaxis_title='x',
yaxis_title='y',
hovermode='x unified',
sliders=sliders,
hoverlabel=dict(bgcolor='white', font_size=12)
)
fig.show()
Выводы¶
Построенные графики и их близость к точному решению полностью согласуются с погрешностями методов, с помощью которых эти графики решений построены: чем точнее метод, тем ближе решение, полученное с его применением к точному решению задачи Коши